Skyrmion versus vortex flux lattices in p-wave superconductors 
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p-wave superconductors allow for topological defects known as skyrmions, in addition to the usual 
vortices that are possible in both s-wave and p-wave materials. In strongly type-II superconductors 
in a magnetic field, a skyrmion flux lattice yields a lower free energy than the Abrikosov flux lattice 
of vortices, and should thus be realized in p-wave superconductors. We analytically calculate the 
energy per skyrmion, which agrees very well with numerical results. From this, we obtain the 
magnetic induction B as a function of the external magnetic field H, and the elastic constants of 
the skyrmion lattice, near the lower critical field Hd- Together with the Lindemann criterion, these 
results suffice to predict the melting curve of the skyrmion lattice. We find a striking difference 
in the melting curves of vortex lattices and skyrmion lattices: while the former is separated at 
all temperatures from the Meissner phase by a vortex liquid phase, the skyrmion lattice phase 
shares a direct boundary with the Meissner phase. That is, skyrmions lattices never melt near Hd, 
while vortex lattices always melt sufficiently close to Hd- This allows for a very simple test for 
the existence of a skyrmion lattice. Possible /^SR experiments to detect skyrmion lattices are also 
discussed. 

PACS numbers: 74.50.+r,74.70.Pq,74.25.Fy 



I. INTRODUCTION 



One of the most fascinating phenomena exhibited by 
conventional, s-wave, type-II superconductors is the ap- 
pearance of an Abrikosov flux lattice of vortices in the 
presence of an external magnetic field H in a, range 
Hci < \H\ < Hc2 between a lower critical field Hd and 
an upper critical field Hc2 j^ It has been known for quite 
some time both theoretically^i^'^'^ and experimentally^ 
that these flux lattices can melt. The melting curve sep- 
arates an Abrikosov vortex lattice phase from a vortex 
liquid phase, and the vortex lattice is found to melt in 
the vicinity of both Hd and Hc2, as shown in Fig. [T] 
The melting occurs because the elastic constants of the 
flux lattice (i.e., the shear, bulk, and tilt moduli) van- 
ish exponentially near these fleld values. As a result, in 
clean superconductors, root-mean-square positional ther- 
mal fluctuations ^/{\u(x)\^ grow exponentially as these 
fields are approached. According to the Lindemann cri- 
terion, when these fluctuations become comparable to 
the lattice constant a, the translational order of the flux 
lattice is destroyed; i.e., the lattice melts. 

Vortices are topological defects in the texture of the su- 
perconducting order parameter, and in s-wave supercon- 
ductors, where the order parameter is a complex scalar, 
only one type of defect is possible. In p-wave supercon- 
ductors, the more complicated structure of the order pa- 
rameter allows for an additional type of topological defect 
known as a skyrmion. In contrast to vortices, skyrmions 
do not involve a singularity at the core of the defect; 
rather, the order parameter field is smooth everywhere, 
as illustrated in Fig. [2 Skyrmions were first introduced 
in a nuclear physics context by Skyrme,^ and slight vari- 
ations of this concepli^ were later shown or proposed to 
be important in superfluid ^He )^°'^^ in the Blue Phases of 
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FIG. 1: External field (H) vs. temperature (T) phase diagram 
for vortex flux lattices. Shown are the Meissner phase, the 
vortex lattice phase, the vortex liquid, and the normal state. 
Notice that the vortex lattice is never stable sufficiently close 
to Hd- 



liquid crystals,^- in Quantmn Hall systemSfi^ii^ in itin- 
erant ferromagnets,^^ and in p-wave superconductors^— 
In the latter case, skyrmions carry a quantized magnetic 
flux, as do vortices, although the lowest energy skyrmion 
contains two flux quanta, while the lowest energy vor- 
tex contains just one. For strongly type-II superconduc- 
tors, skyrmions have a lower free energy than vortices, 
and a vortex lattice should thus be the state that occurs 
naturallyji^ 

Recent evidence of p-wave superconductivity in 
Sr2Ru0 4^^i^^i^^ provides a motivation for further ex- 




determine the elastic properties of the skyrmion lattice, 
and we predict the magnetic induction distribution n(B) 
as observed in a /iSR experiment. 



II. FORMULATION OF THE SKYRMION 
PROBLEM 



FIG. 2: Order parameter configurations showing a vortex (a), 
and a skyrmion (b). The local order parameters are repre- 
sented by arrows on loci of equal distance from the center of 
the defect. If the order parameter space is two-dimensional, 
only vortices are possible, and there is a singularity at the 
center of each vortex, (a). If the order-parameter space is 
three-dimensional, a skyrmion can form instead, where the 
spin direction changes smoothly from "down" at the center 
to "up" at infinity, (b). 



In this section we review the formulation of the 
skyrmion problem presented in Ref. Ilg . who derived an 
effective action that allows for skyrmions as saddle-point 
solutions. The resulting ordinary differential equations 
(ODEs) describing skyrmions^^ are the starting point for 
our analytic treatment. 



A. The action in the London approximation 



ploring the properties of skyrmion flux lattices in such 
systems.^- It was shown numerically by Knigavko et al. — 
that the interaction between skyrmions falls off only as 
1/R with distance R, as opposed to the exponentially de- 
caying interaction between vortices. As result, skyrmion 
lattices have a very different dependence of the magnetic 
induction on the external magnetic field near Hd than 
do vortex lattices. In this paper we confirm and ex- 
pand on these results. We show analytically that the 
skyrmion-skyrmion interaction, in addition to a lead- 
ing l/i?-dependence, has a correction proportional to 
XiiR/B? that explains a small discrepancy between the 
numerical results in Ref.KLg and a strict \/R fit, and we 
calculate the interaction energy up to 0{1/R^). We fur- 
ther show that the melting curve of a skyrmion lattice 
is qualitatively different from that of a vortex lattice. 
Namely, skyrmion lattices melt nowhere in the vicinity 
of TJci, so there is a direct transition from the Meiss- 
ner phase to the skyrmion lattice, see Fig. [5] below. Fi- 
nally, we predict and discuss the magnetic induction dis- 
tribution n{B) of a skyrmion lattice state as observed in 
a muon spin resonance (/xSR) experiment. For a vor- 
tex lattice, the exponential decay of the magnetic in- 
duction B at large distances from a vortex core implies 
n{B) oc XuB/B. For a skyrmion lattice, we find that B 
decays only algebraically, which leads to n{B) oc B^^l^ . 
Some of these results have been reported before in Ref. 

The paper is organized as follows. In Sec. [IT] we re- 
view the formulation in Ref. 'lo of the skyrmion prob- 
lem. In particular, we start from the Ginzburg-Landau 
(GL) model for p-wave superconductors and consider the 
free energy in a London approximation. We parameterize 
the skyrmion solution of the saddle-point equations, and 
express the energy in terms of the solution of the saddle- 
point equations. In Sec. IIIII we analytically solve these 
saddle-point equations perturbatively for large skyrmion 
radius i?, and we calculate the energy of a single skyrmion 
as a power series in 1/i? to order \jR?. In Sec. IIVI we 



We start from a Landau-Ginzburg-Wilson (LGW) 
functional appropriate for describing spin-triplet super- 
conducting order. 



dx L{ii{x),A{x)), 



(2.1a) 



with an action density 



-^(^^^)^- 



an 

(2.1b) 



Here t\)(x) is a 3-component complex order parameter 
field, ^^ A(x) is the electromagnetic vector potential, and 
Z> = V — iqA denotes the gauge invariant gradient oper- 
ator, m and q are the mass and the charge, respectively, 
of a Cooper pair, and we use units such that h = c = \. 
t, u, and V are the parameters of the LGW theory. 

Let us look for saddle-point solutions to this action. 
In a large part of parameter space, namely, for v < 
and u > —V, the stable saddle-point solution is has the 
form ip{x) = tp = /o (l,i,0)/-\/2, where the amplitude 
/o is determined by minimization of the free energyi^ 
This is known as the /3-phase, and it is considered the 
most likely case to be realized in any of the candidates 
for p-wave superconductivity^. Fluctuations about this 
saddle point are conveniently parameterized by writing 
the order parameter field as 



ipix) 
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f{x){n{x) + im{x)), 



(2.2) 



where n{x) and m{x) arc unit real orthogonal vectors 
in order-parameter space and f{x) is the modulus of or- 
der parameter. With this parameterization, the action 
density can be written 
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(n ■ difh 



(V X A)2, 



(2.3) 



where I = n x rh, summation over repeated indices is 
implied, and we have made use of the identities Usted in 
Appendix [A] 

There are two length scales associated with the action 
density, Eq. (|2.3|) . The coherence length ^ is determined 
by comparing the /^ term with the (V/)^ term. 



e = l/V2m|i 



(2.4a) 



It is the length scale over which the amplitude of the 
order parameter will typically vary. The London pene- 
tration depth A is determined by comparing the A^ term 
with the (V x A)'^ term, 



A= V"^/47rq2(/)5 



(2.4b) 



The ratio of these two length scales, k = \/^, is the 
Ginzburg-Landau parameter. Now we write f{x) = 
/o + Sf{x), with /o — ^/—t/2{u + v). Deep inside the 
superconducting phase, where — t > is large, the ampli- 
tude fluctuations Sf are massive, and to study low-energy 
excitations one can integrate out / in a tree approxima- 
tion. This approximation becomes exact in the limit of 
large n and is known in this context as the London ap- 
proximation. We introduce dimensionless quantities by 
measuring distances in units of A and the action in units 
of $Q/327r'^A, and we introduce a dimensionless vector 
potential a = 2ttXA/^q, with $o = 27r/g the magnetic 
flux quantum. Ignoring constant contributions to the 
action we can then write the action density in London 
approximation as follows, ^^ 



cl = 1 {d.rif 



{fiditfi - aif + b^ 



(2.5) 



with b = V X a. The above derivation makes it clear that 
this effective action is a generalization of the 0(3) non- 
linear sigma model (represented by the first term on the 
right-hand side of Eq. (|2.5p ) that one obtains for a real 3- 
vector order parameter by integrating out the amplitude 
fluctuations in tree approximationi^^ 



B. Saddle-point solutions of the effective action 

We now are looking for saddle-point solutions to the 
effective field theory, Eq. (|2.5|) . Considering I and n in- 
dependent variables, and minimizing with respect to I 
subject to the constraints l^ — it? — 1 and Z • n = 
yields 



with 



VH - l{l ■ V^Z) + 2Ji{l X d,l) = 0, (2.6a) 



J = V X 6 (2.6b) 



the supercurrent. The variation with respect to a is 
straightforward and yields a generalized London equa- 
tion. 



It is convenient to take the curl of Eq. (|2.6cP and use 
Eq. (|A3p to express the right-hand side of the resulting 
equation in terms of I. We then obtain the saddle-point 
equations as a set of coupled partial differential equations 
(PDEs) in terms of b and I only: 



2 1 



(djl X dkl), (2.7a) 



V^Z - l{l ■ V^Z) + 2eijk djbk{l x dd) = 0. (2.7b) 



Notice that the right-hand side of Eq. (j2.7ap is valid in 
this form only at points where l{x) is differentiable, see 
Eq. (JA3|. Field configurations that obey these PDEs 
have an energy 



E^ dx 



]:{^^if + {n ■ d^m - af + b^ - 2h ■ b 



(2.8) 

where we have added a uniform external magnetic field 
h measured in units of <I>o/27rA^. Notice that the energy 
depends on n and m, whereas Eqs. ()2.7p depend only on 
Z, and that different choices of n and m can lead to the 
same I. Therefore, a field configuration satisfying Eqs. 
(j2.7p is only necessary for making the energy stationary, 
but not sufficient. 



1. Meissner solution 

A very simple order parameter configuration consists 
of constant h{x) and m{x) everywhere, see Fig. [S] This 




Qi + Ji — hdirh. 



(2.6c) 



FIG. 3: Configurations of the vectors I, rh, and n in a Meiss- 
ner phase. All three vectors point in the same direction ev- 
erywhere. 

leads to an l{x) = I that is constant everywhere. Equa- 
tion (j2.7bp is then trivially satisfied. The right-hand side 
of Eq. (|2.7ap vanishes, and hence the PDE for b reduces 
to the usual London equation with a solution b(x) = in 
the bulk. This solution describes a Meissner phase with 
energy £;m = 0. 



2. Vortex solution 

Now consider a field configuration where ^{x) and 
rh^x) are confined to a plane (say, the x-y plane), but 
rotate about an arbitrarily chosen point of origin: 



n(a;) = {cos(j),sm(j),0), 
7ri{x) = (— sini/), COS0, 0), 



(2.9) 



where cf) denotes the azimuthal angle in the x-y plane with 
respect to the x-axis. This field configuration, known as 
a vortex and shown in Fig. 31 corresponds to a constant 
I everywhere except at the origin, where there is a singu- 
larity. Therefore, the right-hand side of Eq. (|2.7aP is not 




FIG. 4: Configurations of the vectors I, rh, and h for a vortex. 
£ is constant, whereas m and n rotate about the vortex core. 
Notice that the vector shown in Fig. [2ja) is h. 



applicable, and we return to Eq. ()2.6cp . which takes the 
form 



flj + ^iikdjUk = di> 



(2.10) 



For any closed path C in the x-y plane that surrounds 
the origin one has 

<j> de.-V(t){x) = 2iT, (2.IIa) 

or, by Stokes' theorem, 

/ ds • (V X V4>{x)) = 27r, (2.11b) 

J A 

where A is the surface whose boundary is C?^ This quan- 
tization condition shows that, instead of Eq. (|2.7ap . we 
have 



bix) - V'^b{x) = 27rz(5(x) Sijj). 



(2.12) 



This is solved by a b that is equal to the boundary con- 
dition value everywhere along the z-axis and that falls 
off exponentially away from the z-axis. This solution is 
known as a vortex, and the amount of magnetic flux con- 
tained in one vortex is one flux quantum <i>o-^^ It is the 



precise analog of, and, indeed, essentially identical to, the 
familiar vortex in conventional s-wave superconductors. 

The energy of a vortex given by Eq. (|2.12p . as calcu- 
lated from Eq. (|2.8p . is logarithmically infinite. This is 
due to the point-like nature of the vortex core where the 
amplitude of the order parameter goes discontinuously 
to zero. In reality, the amplitude cannot vary on length 
scales shorter then the coherence length ^, which provides 
an ultraviolet cutoff. The energy is then proportional to 
Inn} In an external magnetic field this energy cost is 
offset by the magnetic energy gain due to letting some 
flux penetrate the sample. For n larger than a critical 
value Kc = l/v^, and for external fields larger than the 
lower critical field Hd, a hexagonal lattice of vortices 
has a lower energy than the Meissner phase. This state 
is known as an Abrikosov flux lattice, 

and is precisely the same as that in conventional s-wave 
superconductors li 



3. Skyrmion solution 

Due to the three-component nature of the order pa- 
rameter, more complicated solutions of the saddle-point 
equations can be constructed for which the vector I is 
not fixed. Let 6 be the angle between I and the z-axis, 
and consider a cylindrically symmetric field configuration 
parameterized as 

I — ezCosd{r) -\- e.rSm9{r), 
h = (e^ sin 0(r) — Cr cos 6'(r)) sin iy9 -I- e^ cos C(9 
rh = {ezSnv6{r) — €.rCos9{r)) costf ~ e.^pSV[\^p. 

(2.13) 

For this to minimize the energy, I at large distances from 
the origin must be constant because of the first term in 
the energy, Eq. (|2.8p . and for a skyrmion centered in a 
cylinder of radius R we take I to point in the +z-dircction 
for r — R, d{r = R) = 0. The quantization condition 
analogous to Eq. (j2.11bp for the vortex isS^^ 



dxdy Eij I ■ (Oil X djl) — Stt 



(2.14) 



To be consistent with this, I must point in the —z- 
direction at the origin, 0{r — 0) = n. 

Equation (|2.13p parameterizes the order parameter in 
terms of a function 9{r). In addition, the energy depends 
on the vector potential which we take to be purely az- 
imuthal, in accordance with our cylindrically symmetric 
ansatz, 



a{x) = a{r) 



(2.15) 



With this parameterization, we obtain from Eq. ,, ,, 

the energy per unit length, along the cylinder axis, of a 



cylindrically symmetric skyrmion in a region of radius R, 



E/Eo 



drr 



dr 'I 



dr 'I 



{9'{r)y 



sin^ 9{r) 



1 2 



- (1 + cos 6'(r)) +a{r) 
r 



a{r) 



a'{r) 



(2.16) 



where Eq = (J>o/47rA)^. This expression was first ob- 
tained in Ref. ^^ The three terms correspond to the 
three terms in the London action, Eq. ()2.5|) . They repre- 
sent the energy of the nonhnear sigma model, the kinetic 
energy of the supercurrent, and the magnetic energy, re- 
spectively. Minimization of E with respect to 9{r) and 
a{r) yields Euler-Lagrange equations 



r{r) + -e'{r) 



sin6'(r) 



2 + cos 6l(r) 



2a{r) 



a"{r) + - a'{r) ^ "('') = «(^) + - [1 + cos 6l(r)] . 



(2.17a) 

r)]. 

(2.17b) 



This set of coupled, nonlinear ODEs must be solved 
subject to the boundary conditions 9{r = 0) = tt and 
9{r = R) = 0, as explained above. The solution is 
known as a skyrmion, and each skyrmion contains two 
flux quantaj^ Since Eqs. (|2.7p are necessary for making 
the energy stationary, the solution of Eqs. (|2.17p . inserted 
in Eqs. (|2.13l[^TT5)) . is guaranteed to be a solution of Eqs. 
((2Jl) as well. 

The energy of a single skyrmion is finite even in London 
approximation, see Sec. lIIII below. For large values of the 
Ginzburg-Landau parameter k a skyrmion therefore has 
a lower energy than a vortex, and the value of the lower 
critical field Hd, at which the Meissner phase becomes 
unstable, is correspondingly lower for skyrmions than for 
vortices. This is the basis for the expectation that, in 
strongly type-II (i.e., large-K) p-wave superconductors, a 
skyrmion flux lattice will be realized rather than a vortex 
flux lattice. 



III. ANALYTIC SOLUTION OF THE 
SINGLE-SKYRMION PROBLEM 



use Eq. (j2.17bp to iteratively express a in terms of 9 and 
its derivatives. Substitution in Eq. (|2.17ap then yields a 
closed ODE for 9{r) that has to be solved. 



Zeroth order solution 



Let us first consider R — oo. For 



cx), the left-hand 



side of Eq. (|2.17bp falls off as l/r^, and hence the vector 
potential, to zeroth order for large r, is given by 



aoc{r) == -- [1 + cos 6'(r)] 
r 



(3.1) 



Note that we use the exact ^(j') in this expression, not 
the zeroth order approximation to it. Since we can only 
compute 9{r) perturbatively, this expression for the ze- 
roth order vector potential will itself have to be expanded 
perturbatively later. Substitution in Eq. (|2.17aD yields 

r^ 9"{r) +r9'{r) = ]- sm{29{r)). (3.2) 

The solution obeying the appropriate boundary condition 

;„16 



with 



.{r)^I{r/(.), 



f[x) — 2arctan(l/a;). 



(3.3a) 



(3.3b) 



The length scale € is arbitrary at this point and will be 
determined later from the requirement 9{r = R < cx)) = 
0. For _R ^ 1 it will turn out that (. ex \/i?- The skyrmion 
solution is schematically shown in Fig. [5l 




We now need to solve the coupled ODEs (|2.17p . Due 
to their nonlinear nature, this is a difficult task, and in 
Ref. Il6| it was done numerically. It turns out, however, 
that one can construct a perturbative analytical solution 
in the limit of large skyrmion radius, i? S> A, with X/ R 
as a small parameter. This provides information about 
the superconducting state near i?ci, where the system is 
always in that limit. We will construct the perturbative 
solution, and calculate the energy, to second order in the 
small parameter. Our general strategy is as follows. We 



FIG. 5: Configurations of the vectors I, m, and n for a 
skyrmion. Notice that the vector shown in Fig. [IJb) is l. 



B. Perturbation theory for R^ \ 

We now determine the corrections to the zeroth or- 
der solution. Let us write 9{r) = 9oo{r) + d9{r) and 



a(r) = aoo(r') + Sa{r) and require |(5a(r)| ^ |aoo('')| and 
\d9{r)\ < 1^ An inspection of the ODEs (PTT]) shows 
that for r < 0{£) , the corrections can be expanded in a 
series in powers of l/£, 



S9{r) 
Sa{r) 



1 



girl I) - 
a{r/t} 



1 

¥ 
1 

£5 



h{r/t) - 
p{r/i) 



0{\/t^), (3.4a) 
^0{\/f). (3.4b) 



The functions a and (3 can be determined by substitut- 
ing Eq. (|3.4bp in Eq. (|2.17bp and equating coefficients of 
powers of l/£. The resuhing equations for a and (3 are 
hnear algebraic equations, not ODE's, because terms in- 
volving derivatives of a and (3 only enter at higher order 
in 1/^, as one can verify by direct calculation. Hence, the 
solutions for a and (3 can be read off at once, and are: 



(3{x) 



162: 



= 2 



(l + a;2)3' 
(3a;4 - 6x2 



1) 



9{x) 



x2(l+2;2)3 
-^g"{x)+A{x), 

l_ ry^ii. 



(3^2 - 1) 

a;(l + a;2)2 



(3.5a) 



(3.5b) 



where 



A(x) = a" [x) H — a (x) 5- a{x) 

X X 

384x(a;2 - 1) 



(3.5c) 



Similarly, by comparing coefhcients in Eq. (j2.17ap we find 
ODEs for the functions g and /i, 



g"[x)^-g\x)-\co^{2!{x))g{x) = 
X x^ 

h"{x) + - h'{x) - \ cos{2f{x)) h{x) = 



X 



x^ 



i{f{x))a{x) 



(3.6a) 



sin(/(x))/3(x) J sin(2/(a;)).g (x) cos{f{x))a{x)g{x), 



(3.6b) 



r 



with f{x) from Eq. ((3^ . 

The ODE (|3.6aP for g can be solved by standard meth- 
ods, see Appendix |B1 The physical solution is the one 
that vanishes for a; ^ 0; it is proportional to x for x S> 1. 
We find 



lix) 



4 x[x^{4 + x2) + 2(1 + x^) ln(l + x'^)] 

3 (1+X2)2 



(3.7a) 



the large-x asymptotic behavior of which is 



4 16 Inx 8 

gix»l) = --x-- — -- 



O(i^). (3.7b) 



This determines both the function P{x), Eq. (j3.5bp . 
and the inhomogeneity of the ODE (j3.6bp for h{x). The 
latter can again be solved in terms of tabulated functions, 
see Appendix [Bl but we will need only the two leading 
terms for x — *■ 00. The physical solution is again the 
one that vanishes for x —> 0, and its large- a; asymptotic 
behavior is 



h{x > 1) 



32 



xlnx ■ 



536 
135 



x + 0{l/x). (3.8) 



Finally, we need to fix the length scale £. It is deter- 
mined by the requirement 9{r = R) ^ 0. We find 



t = ^^R 



- — In R:^ 

R R 



O 



flnR 

VrN^ 



where 



and 



12 



(7-6d/c2)-- ln(c/2) 



c = 4/3, 

d = 536/135, 



,(3.9a) 



(3.9b) 



(3.9c) 
(3.9d) 



are the absolute values of the coefhcients of the terms pro- 
portional to X in the large- a; expansions oi g{x) and h{x), 
respectively. We see that, for i? ^ 1, ^ is indeed propor- 
tional to \fR , as we had anticipated above. That is, the 
characteristic skyrmion length scale i is the geometric 
mean of the London penetration depth A (recall that we 
measure all lengths in units of A) and the skyrmion size 
R. We now can also check our requirement 59 <^\: from 



Eq. (|Oa|) we see that for r < £, 59 {r) ex 1/R, while 
for r ^ £, 59{r) is bounded by a term proportional to 
1/i?^/^. For R large compared to the penetration depth 
the condition is thus fulfilled for all r. Similarly, 5a is 
found to be small compared to Ooo for all r. 



C. Energy of a single skyrmion 



By using our perturbative solution in Eq. (|2.16p . we 
are now in a position to calculate the energy of a single 
skyrmion to 0{1/R^). It is convenient to first expand 
the energy in powers of 1/^^, and then determine the 
i?-dependence by using Eqs. (|3.9p . 

Let us first consider the supercurrent energy Ec, i.e., 
the second term in Eq. (|2.16p . It can be written 

.R -, .R 

Ec/Eo=l drr{5a{r)f = j^l drr {a{r/i)f 

+ 0(1/ f). (3.10) 

Using Eqs. (|3.5ap we find 



39 1 

E,/E,^-j^ + 0{l/i'). 



(3.11) 



Now consider the magnetic energy E^, which is the 
third term in Eq. (|2.16p . It can be written 



Em/Eo 



drrb (r). 



(3.12) 



with 



h{r) = - a^{r) + a'^{r) + - 5a{r) + 5a'{r), (3.13a) 
r r 

the magnetic induction in our reduced units. Notice that 
in calculating aoo{r)^ 9{r) in Eq. p.ip needs to be ex- 
panded to first order in 59^ as noted earlier. The two 
leading contributions to b^ are then 

^2 



&nO - ^ 



1 



1 



16 

¥ (l + a;2)4 "^ (l + 2;2) 

^ i^ +-a{x)+a'{x) 



1 



■9{x) 



x{\ + x'^ 
0(l/£^), (3.13b) 



(3.14) 



1 + x^ X 
where x = r/L Performing the integral yields 

Finally, we need to calculate the energy Eg coming 
from the gradient terms in the first term in Eq. (|2.16p . 
The expansion of the two terms in the integrand yields 
seven integrals that contribute to the desired order, they 
are listed in Appendix [Cl The result is 
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FIG. 6; Numerical data for the energy per skyrmion per unit 
length (circles) together with the best fit to a pure 1/R behav- 
ior (dashed line) from Ref.^^, and the perturbative analytic 
solution given by Eq. p.l6p (solid line). A numerical so- 
lution using spectral methods is indistinguishable from the 
perturbative one. 



Adding the three contributions, and using Eqs. p.9p . 
we find our final result for the energy of a skyrmion of 
radius i? :^ 1, 



E/Eo = 2 



8\/6 1 
~3~^ 



16 InR 
y "R2- 



-^[7 + 301n(3/2)]ij 



0{ln^R/R^). 

(3.16) 



Knigavko et alf^^ solved the Eqs. (|2.17p numerically, 
and thereby numerically determined the energy, which 
they fit to a l/i?-dependence. Their results are shown 
in Fig. [6] together with the analytical result given in Eq. 
p.l6p . The perturbative solution up to 0(\nR/R'^) was 
first given in Ref. [2l|. We have also solved the equa- 
tions numerically, using spectral methods to convert the 
boundary value problem to a set of algebraic equations 
for the unknown coefficients in an expansion in Cheby- 
shev polynomials^^. For the i?-range shown, and on the 
scale of the figure, the result is indistinguishable from the 
perturbative one. 



OBSERVABLE CONSEQUENCES OF THE 
SKYRMION ENERGY 



IV. 



Our calculation of the skyrmion energy in Sec. Illll has 
been for a cylindrically symmetric skyrmion. The result 
shows that each skyrmion will try to maximize its ra- 
dius in order to minimize the energy, which leads to a 
repulsive interaction between skyrmions whose potential 
is proportional to 1/R. Skyrmions are thus expected to 
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form a lattice structure, as do vortices, and they will 
thus not be cylindrically symmetric, since the lattice is 
not. One expects a hexagonal lattice, as in the case of the 
vortex lattice, and our treatment involves the same ap- 
proximation as in the numerical work of Ref. |l6|; namely, 
approximating the hexagonal unit cell by a circle of the 
same area. We expect this approximation to recover the 
correct scaling of the energy, and to reproduce the coef- 
ficients of that scaling to the same accuracy as radius of 
the circle of the same area reproduces the distance from 
the center of a hexagon to the nearest point on its edge; 



i.e., \j2^/3/n — 1 w 0.05. We will now proceed to cal- 
culate observable consequences of the dependence of the 
energy on the radius of the unit cell. These include the 
relation B{H) between the magnetic induction B and 
the external magnetic field H , the elastic properties of 
the skyrmion lattice and the resulting phase diagram in 
the _ff-T-plane, and the /iSR signature of the skyrmion 
lattice. 



A. B(H) for a skyrmion lattice 

We start by calculating the dependence of the equilib- 
rium lattice constant R on an external magnetic field H. 
This is done by minimizing the energy per unit volume, 
which is the energy per unit length per skyrmion, Eq. 
(|3.16p . divided by the area per skyrmion, ttR^, plus a 
reduction in the energy of ~2^oH/4:Tt due to the exter- 
nal field. The latter is obtained from the last term in Eq. 
(|2.8p by noting that the magnetic flux J dxdy (z-b) = 2$o 
for each skyrmion in the lattice. This negative external 
field contribution must also be divided by ttR^ to give the 
energy per unit volume. Returning to ordinary units, we 
thus find a Gibbs free energy per unit volume 






A 
R^ 



4\/6A 



1^^^) 



where K = $o/27rA^, and 



A=l -H/H,i, 



(4.1) 



(4.2) 



with Hci = K/2^Q. For H < Hd , we have A > 0, and 
the free energy is minimized by i? — oo; i.e., the skyrmion 
density is zero. This is the Meissner phase. For H > Hd 
the free energy is minimized by 



R = Ra = 2\/6A/A, 



(4.3) 



and there is a nonzero skyrmion density. We see that Hd 
is indeed the lower critical field. Note that the equilib- 
rium flux lattice constant Rq diverges as 1/A, whereas 
in the case of a vortex lattice it diverges only logarithmi- 
cally as ln(l/A).^ For the averaged magnetic induction 
B = 2<i>o/7ri?o this imphes 



BiH) 



1 



■ Hd A^ 



(4.4) 



For H -^ Hd from above, B{H) in the case of a skyrmion 
lattice thus vanishes with zero slope, whereas in the case 
of a vortex lattice it vanishes with an infinite slope.— This 
result, with a slightly different prefactor, was first ob- 
tained from the aforementioned numerical determination 
of E{R) in Ref. 16. Note that the only material param- 
eter that appears in this expression for B is -ffci- 



B. Elastic properties of the skyrmion lattice 

Now we turn to the elastic properties of skyrmion lat- 
tice. Let the equilibrium position of the i^^ skyrmion 
line be described by a two-dimensional lattice vector 
Ri = {Xi, Yi), and the actual position by 

ri{z) = {X^ + Ux{Ri, z), Yi + Uy{Ri,z),z), (4.5) 

where u — {ux,Uy) is the two-dimensional displacement 
vector, and we use z as the parameter of the skyrmion 
line. The strain tensor Uq,^ is defined as 



Ua/six) 



1 f dua dup 



2 \ dxp dxa 



(4.6) 



For a hexagonal lattice of lines parallel to the z-axis, the 
elastic Hamiltonian is^^ 



Hal = - I dx [2fl {Ual3{x)Uap{x)) + Xi, {Uaa{x)Y 



+Ktiit\dMx)f 



(4.7) 



Here summation over repeated indices is implied, /i, Al, 
and Ktm are the shear, bulk, and tilt moduli, respec- 
tively, of the lattice, and we now need to determined 
these elastic constants. 

The combination jj. + Xj^ can be obtained by consider- 
ing the energy change of the system upon a dilation of 
the lattice. Let R change from Rq to i?o(l + e), with a 
dilation factor e ^ 1. Such a dilation corresponds to a 
displacement field u{x) = ex±, where x± is the projec- 
tion of X perpendicular to the z-axis.'^"'^ The strain tensor 
is thus Uq/3 = eSap. Inserting this in the elastic Hamil- 
tonian, Eq. (|4.7p . yields the energy per unit volume for 
the dilation, 



Eda/V^2{^i + XL)e^ 



(4.8a) 



This should be compared with the energy as given by Eq. 

63), 



Edii/V = g{Roil + e))- 



KA^ 



5(^0) = 7; \ln^ 



2 KdR^J 



{eRo? 



Ro 



(4.8b) 



967r2A2 ■ 
Comparing Eqs. (I4.8a|) and (j4.8bp yields 

lj. + Xi^^K A^/ 19277^ X^. (4.9) 



To obtain fj, (or Al) separately, we should consider 
shear deformations, which change the shape, but not the 
area, of the unit cell. Since we have already approxi- 
mated the hexagonal unit cell by a circle, this is difficult 
to do, and we resort to the following heuristic method, 
which will give the correct scaling of /i with A (but not 
the correct prefactors). To this end we observe that our 
result for the Gibbs free energy, Eq. (|4.ip . is of the same 
form we would have obtained if the skyrmions interacted 
via a pair potential U{r) that for distances r < Rq is of 
order KX/r, and for larger distances falls off sufficiently 
rapidly that only nearest-neighbor interactions need to 
be considered. Treating the skyrmion lattice as if such 
an "equivalent potential" were the origin of the skyrmion 
energy allows us to calculate the shear modulus as fol- 
lows: 

If the lattice is subjected to a uniform x — y shear - i.e., 
a displacement field u{x) — 2eyx - for which u^y — Uyx — 
e, and all other components oiuap — 0, the elastic energy, 
Eq. (14. 7p predicts an elastic energy per unit volume of 



E/V 



2^6^. 



(4.10) 



Such a shear skews each fundamental triangle of the 
skyrmion lattice by displacing the top (or bottom, for 
the downward-pointing triangles) to the right (or the left, 
for downward-pointing triangles) by an amount of order 
ei?o, where i?o is the skyrmion lattice spacing found ear- 
lier, Eq. (li?^ (see Fig. [T]). This shortens the length of 




FIG. 7: Shearing of the skyrmion lattice results in a change 
in the distance between skyrmion centers, and hence in their 
effective interaction. See the text for additional information. 

one bond of the triangle by an amount of order ei?o , and 
increases the opposite bond's length by the same amount. 
Hence, the linear in e change in the "equivalent poten- 
tials" of these two bonds cancels, and the total change 
(Ai?/triangle) in the energy per unit length of funda- 
mental triangle, per triangle, is given by: 



AS 
triangle 



U"{Ro){eR^f X 0(1) 



(4.11) 



where the 0(1) factor includes both geometrical factors 
(e.g., sines and cosines), and counting factors (e.g., to 
avoid multiple counting of each triangle). If we take 



U{r) ~ K\/r as suggested above, we have 

K\ 
U"{Ro) = -^ X 0(1). 



Inserting this into Eq. (|4.1ip gives 
AE K\ 



triangle i?o 



e' X 0(1). 



(4.12) 



(4.13) 



This is the change in energy per unit cell. To get the 
energy per unit volume, we must divide by the unit cell 
area, which is ttR^. Doing so gives 

AE _ K\ 

Comparing this with Eq. (j4.10p then determines /i 
K\ 



(4.14) 



M 



Rl ^ ^(^)- 



(4.15) 



Using Eq. (|4.3p for Rq then leads to our final result for 
/i: 



A* 



A2 



xO(l) 



(4.16a) 



From Eq. (|4.9p we see that the bulk modulus or Lame 
coefficient is given given by the same expression. 



At 



^xOd), 



(4.16b) 



We now turn to the tilt modulus i^tiit- This can 
be obtained by considering a uniform tilt of the axes 
of the skyrmions away from the z-axis, i.e., away from 
the direction of the external magnetic field iJ, by an 
angle d <^ 1. For small ■&, "d — \du/dz\. Therefore, 
the tilt energy in Eq. (|4.7p is identical with the change 
of the B ■ H term in Eq. (|2.8p . This contribution to 
the energy is, per unit length and in ordinary units, 
given by —^oHcos9/2tt, and its change due to tilting 
is $0^(1 - cos ^)/27r« <^oHd^/4:Tr = ^aH \d^u\^/ATT. 
Dividing this result by the unit cell area ttRq, using Eq. 
(|4.3p for Rq , and identifying the result with the tilt term 
in the elastic Hamiltonian, Eq. (|4.7p . yields Ktm in the 
vicinity of Hd , 



iZTT 



(4.17) 



We now are in a position to calculate the mean-square 
positional fluctuations (|tt(a;)p). Taking the Fourier 
transform of Eq. (|4.7p , and using the equipartition theo- 
rem, yields 



H^)\')t 



kv.T 



X^ 



^ ,t^z/^'?l + ^"i"?' 



for the transverse fluctuations, and 
kY^T V- 1 



i«(^)r>L = 



V ^^^i2^i + \L)ql+Kultq^, 



(4.18a) 



(4.18b) 
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for the longitudinal ones. Here q± and g^ are the projec- 
tions of the wave vector q orthogonal to and along the 
2-direction, respectively. The Brillouin zone BZ of the 
skyrmion lattice is a hexagon (which we have approxi- 
mated by a circle) of edge length O(l)/i?o in the plane 
perpendicular to the z-axis, and extends infinitely in the 
z-direction. 

Since /i and Al are the same apart from a prefactor of 
0(1) which we have not determined, see Eqs. (|4.16p . the 
same is true for the transverse and longitudinal contri- 
butions to the fluctuations, and it suffices to consider the 
former. Performing the integral over qz yields 



l«(^)l' 






(l«(^)p> 

fceT 
VM-^tiit 



BZ 



87r2 



(4.19) 



The remaining integral over the perpendicular part of the 
Brillouin zone is proportional to I/Rq, and using Eqs. 
(|il^ and (1131) we obtain 



\u{x) 



k^T 1 



XH^^ A3/2 



0(1). 



(4.20) 



Using Eq. (|4.3p again we see that, near Hci, {\u{x)\'^) ex 
Rq <C i?o- That is, in this regime the positional fluctu- 
ations are small compared to the lattice constant, which 
tells us that the lattice will be stable against melting. To 
elaborate on this, let us consider the Lindemann crite- 
rion for melting, which states that the lattice will melt 
when the ratio Fl = {\u{x)\'^) / R^ exceeds a critical value 
Fc — 0(1). In our case. 



fcsT 
^c\A5/2 



vl/2 



xO(l). 



(4.21) 



As iJ — > -ffcij A — > 0, and the Lindemann ratio van- 
ishes. Hence, the skyrmion lattice does not melt at any 
temperature for H close to Hd- 

We finally determine the shape of the melting curve 
Hin(T) near the superconducting transition temperature 
Tc- Since, in mean field theory, Hd oc (Tc — T), and 
A oc l/\/Tc — T^ we find from Eq. 
Fl = const. = 0(1), 

H^ - Hd (X (Te - T) 



^^nn\ by putting 



5/2 



(4.22) 



The resulting phase diagram is shown schematically in 
Fig. ISl Comparing with Fig. H] we see the qualitative 
difference between the vortex and skyrmion flux lattices: 
whereas the vortex lattice always melts near 7?ci, the 
skyrmion lattice melts nowhere near Hd- This is a di- 
rect consequence of the long-ranged interaction between 
skyrmions, as opposed to the screened Coulomb interac- 
tion between vortices. 
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FIG. 8: External field (-ff ) vs. temperature (T) phase diagram 
for skyrmion flux lattices. In contrast to the vortex case, see 
Fig. IT] there is a direct transition from the skyrmion flux 
lattice to the Meissner phase. The theory predicts the shape 
of the melting curve only close to Tc, see Eq. (|4.22p : the rest 
of the curve is an educated guess. 



C. /iSR signature of a skyrmion flux lattice 

Muon spin rotation [fiSR] is a powerful tool which 
has been extensively applied to study the vortex state in 
type-II superconductors j 32i33 ^ crucial quantity in this 
type of experiment is the ^SR line shape n{B), which 
is the probability density that a muon experiences a lo- 
cal magnetic induction B and precesses at the Larmor 
frequency that corresponds to B. It is defined as 



n{B)^{d{B{x)-B)), 



(4.23) 



where B(x) is the magnitude of the local magnetic in- 
duction, and (. . .) denotes the spatial average over a flux 
lattice unit cell. 

To predict the ^SR line shape for a skyrmion flux lat- 
tice near Hd it is sufficient, for large Rq, to use only 
the lowest solution for the magnetic induction obtained 
in Sec, nil Al Inserting Eqs. dSSD into Eq. (IXT]) . we find 
for the magnetic induction in reduced units 



6(r) = -- 



4^2 



(r2+£2)2- 
Restoring physical units then gives 
i^ciA^ e 



B{r) 



(r2+<'2)2' 



(4.24) 



(4.25) 



where we've dropped the minus sign since only the mag- 
nitude of B can be detected in /iSR measurements. 

From Eq. (|4.23p we then find, for H near iJ^i, where 
our theory is valid. 



1 /ifciA 



B 



(skyrmions). (4.26) 

Hd 
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Of course, n{B) is only non-zero for those values of B 
that actually occur inside the unit cell of the skyrmion 
lattice. From Eq. (j4.25p . we see that the maximum value 
of B will occur at the center of the unit cell (r = 0), 
which gives 

\BU. = \B{r = 0)\ = ^ = ^, (4.27a) 

The minimum value of B occurs at the edge of the unit 
cell (i.e., r = R), where Eq. (|4.25p gives 



\B\ 



\Bir^R)\ 



H.iXH^ iJelA3 



2i?4 



288 



(4.27b) 



In the second equalities in Eqs. (|4.27p we have used Eqs. 
()3.9|) and ()4.3p to express £ in terms of R and R in terms 
of A, respectively. 

To summarize: the prediction of our cylindrical ap- 
proximation for n{B) is that the simple power law Eq. 
(|i:^ holds for Bmin < B < B^^^. For B < B^i^ or 
B > B„,ax, n{B) = 0. 

Since the above results were derived in the cylindrical 
approximation, we expect the numerical coefficients in 
Eqs. (|4.27|) to be off by the approximately 5% mentioned 
in the opening paragraph of Sec. IIVI throughout most of 
the range Bmin < B < B^^^. When B gets close to 
Bmin, however, we expect more radical departures from 
the cylindrical approximation. This is because contours 
of constant B near the edge of the hexagonal unit cell will, 
for B within 5% or so of Smin or so, start intersecting the 
unit cell boundary, leading to van Hove-like singularities 
in n(B). Such subtleties cannot be captured within the 
cylindrical approximation. Note, however, that they only 
occur over a very small range of B; for the remainder of 
the large window i?niin < B < B^^i^^ (which spans three 
decades even for A as big as 0.2), Eq. (|4.26p holds, up 
to the aforementioned 5% numerical error in its overall 
coefficient. 

To compare this result with the corresponding one for 
a vortex flux lattice, we recall that in that case B{r) is 
given by a modified Bessel function which for distances 
r ^ A takes the form 



B{r) (X 



1 



,-rl\ 



For small -B, we then find from Eq. (|4.23p 

ln(l/S) 



n{B) ex 



B 



(vortices) 



(4.28) 



(4.29) 



We see that the /iSR line shape is qualitatively different 
in the two cases, due to the long-range nature of B{r) in 
the skyrmion case versus the exponential decay in the 
vortex case. 



referred to as skyrmions, rather than by ordinary vor- 
tices. For strongly type-II materials in the /3-phase, 
skyrmions are more stable than vortices.— We have pre- 
sented an analytical calculation of the energy of a cylin- 
drically symmetric skyrmion of radius i? up to 0{1/ R^) 
in an expansion in powers of \/ R. This provides excel- 
lent agreement with numerical solutions of the skyrmion 
equations. The interaction between skyrmions is long- 
ranged, falling off only as the inverse distance, in contrast 
to the exponentially decaying interaction between vor- 
tices. As a result, the elastic properties of a skyrmion flux 
lattice are very different from those of a vortex flux lat- 
tice, which leads to qualitatively different melting curves 
for the two systems. The phase diagram thus provides a 
smoking gun for the presence of skyrmions. In addition, 
the /iSR line width for skyrmions is qualitatively different 
from the vortex case. 

We finally mention two limitations of our discus- 
sion. First, we have restricted ourselves to a discussion 
of a particular p-wave ground state, namely, the non- 
unitary state sometimes referred to as the /3-phase. This 
state breaks time-reversal symmetry and the recently re- 
ported absence of experimental evidence for the latter 
in Sr2Ru043i, suggests to also consider other possible p- 
wave states and their topological excitations, in analogy 
to the rich phenomenology in Helium 3.^^ Second, in a 
real crystalline material, crystal-field effects will invali- 
date our isotropic model at very long distances, and cause 
the skyrmion interaction to fall off exponentially. This 
is the same effect that makes, for instance, the isotropic 
Heisenberg model of ferromagnetism inapplicable at very 
long distances and gives the ferromagnetic magnons a 
small mass. It should be emphasized that this is usually 
an extremely weak effect that is also material dependent. 
Once p-wave superconductivity has been firmly estab- 
lished in a particular material, this point needs to be re- 
visited in order to determine the energy scales on which 
the above analysis is valid. 
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CONCLUSION 



APPENDIX A: PROPERTIES OF ORTHOGONAL 
UNIT VECTORS 



In summary, we have considered properties of a flux 
lattice formed by the topological excitations commonly 



Let h and rh be orthogonal real unit vectors, and I 



n X rh. Then the normalization condition n,n,- 



m^mi 



1 and the orthogonality condition hirhi = imply 



hjdifij — rhjdiihj — 0, 
hj diiJij — — rfij di hj . 



With these relations it is straightforward to show that 

difijdifij + dirhjdiThj — 2{njdirrij){nkdi'fhk) + dtljdilj. 

(A2) 

Finally, in regions where l{x) is differentiable the 
Mermin-Ho relation^ holds, 

I ■ {dil X djl) ~ difi ■ djTh — diih ■ djh. (A3) 
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(This is the solution that vanishes as a; ^ 0. The sec- 
ond solution diverges in this limit.) Now write F{x) = 

(Ala) Fy,{x)G{x), and let y{x) = G'{x). Then y is found to 

(Alb) obey the elementary first-order ODE 



y'{x) +p{x)y{x) = q{x)/Fi,{x), 



with 



p{x) = [2F{^{x)+Fi,{x)/x]/Fi,ix). 



The solution is 



(B3a) 



(B3b) 



APPENDIX B: SOLUTIONS OF THE ODEs FOR 

g AND h 



The functions g and h in Sec. IIIIBI both satisfy an 
ODE of the form (see Eqs. dSU) 



F"{x) + i F'{x) - KnT^.P n^) = 9(x), (Bl) 
X x^(i + x'^y 

with an inhomogeneity q given by the right-hand side of 
Eq. p.6ap or (j3.6b|) . respectively. It is easy to check that 
the corresponding homogeneous equation, obtained from 
Eq. (jBip by putting q{x) = 0, is solved by 



Fi,{x)=x/{\+x^). 



(B2) 



y{x) = e 



— p- fdxp 



C,+ dxqef^^^ 



(B4) 



with Ci an integration constant. A second integration 
yields G{x), and hence F{x) in terms of two integra- 
tion constants. The latter can be determined by re- 
quiring that for small x the solution coincides with the 
asymptotic solution that vanishes as x — > 0. By using 
a power-law ansatz for g and h in Eqs. p.6p we find 
g{x -> 0) = -8a;3 + 0(a;^ and h{x ^ 0) = 256x^ + 0{x'^), 
which suffices to fix the integration constants. For g(x) 
we find the expression given in Eq. (|3.7ap . For h{x) we 
obtain 



h{x) 



1 



{592 + 2x^ [8(-l, 119 + 90a;2 + 286a;^ + 2402;*^ + 30a;®) + 10, 320(1 + x^f] 



(270a;(H-a;2)4) 

+2, 296(-l + x^){l + x^f + 4a;2(i + x^f + 1, 7041na; 



-F321n(l 



-30 



142X-' 

2\ 



276a;'' -^ 171a;*^ ■\- 52a;® - \hx^^ - 15(3 -f x:^){x + x^f ln(l -I- a;^)] 



-l,920a;2(l-f a;2)3Li2(-a;2)} 



r 



(B5) 



with Li the polylogarithm function. The asymptotic be- with 
havior for large x is given by Eq. 



APPENDIX C: CONTRIBUTIONS TO E^ 

By expanding the integrand of the first term in Eq. 
(I2.16p . we can express the energy E^ to OiXjB?) in terms 
of seven integrals, 



^s/^o = ^/,: + 0(l/ 



(CI) 



h = 



= 4 



fl/£ 



dx ■ 



il + x- 



,2^2 ■ 



2 /"^/^ 1 

rfa; ■ 



£2 7g 1 + a;2 \ x^ 



(a;2 - 1) 






(C2a) 



— g{x)-xg'{x)] , 



(C2b) 



(a;4 - 6a;2 + 1) 



x{g'ix)r+ ''^^^^-^-' g'{x) 



(C2c) 



13 



h = 



h = 



h 



2 /-^/^ , 1 



e> 



1 + x^ 



(x^-l) 



— h{x) — xh'{x) 



R/i 



3£^ 



dx 



(xg'{x)h'{x) 



+ ^^^^— ^TT^ g{x)h{x) 







1(1 + 2:2)2 

«/^ (.2 - 1)2 

ax ^ ?rTTr g [X). 



6^8 7o a;(l + a;2)2 



(C2d) 

(C2e) 
(C2f) 
(C2g) 



Evaluating the integrals to 0{l/£^) yields Eq. (|3.15p . 
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